home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / GD-LINF2.LIB < prev    next >
Text File  |  1985-04-03  |  2KB  |  79 lines
  1.  
  2. { -> 220 }
  3. procedure get_data(var t    : ary;        { independedt variable }
  4.            var p    : ary;        { dependent variable }
  5.            var nrow    : integer);    { length of vectors }
  6. var    i    : integer;
  7.  
  8. begin
  9.   nrow:=10;
  10.   for i:=1 to nrow do
  11.     t[i]:=(i+6.0)*100.0;
  12.   p[1]:=1.0E-9;        p[2]:=5.598E-8;
  13.   p[3]:=1.234E-6;    p[4]:=1.507E-5;
  14.   p[5]:=1.138E-4;    p[6]:=6.067E-4;
  15.   p[7]:=2.512E-3;    p[8]:=8.337E-3;
  16.   p[9]:=2.371E-2;    p[10]:=5.875E-2;
  17.   for i:=1 to nrow do
  18.     p[i]:=ln(p[i])    { take log data }
  19. end;        { procedure get_data }
  20.  
  21.  
  22.  
  23. procedure linfit(X,    { independent variable }
  24.          y    : ary;    { dependent variable }
  25.     var y_calc    : ary;    { calculated dep. variable }
  26.     var resid    : ary;    { array of residuals }
  27.     var coef    : arys;    { coefficients }
  28.     var sig        : arys;    { error on coefficients }
  29.         nrow    : integer;    { length of ary }
  30.     var ncol    : integer);    { number of terms }
  31.  
  32. { least-squares fit to nrow sets of x and y pairs of points }
  33. { Seperate  procedure needed:
  34.     SQUARE -> form square coefficient matrix
  35.     GAUSSJ -> Gauus-Jordan elimination }
  36.  
  37. var    xmatr        : ary2;        { data matrix }
  38.     a        : ary2s;    { coefficient matrix }
  39.     g        : arys;        { constant vector }
  40.     error        : boolean;
  41.     i,j,nm        : integer;
  42.     xi,yi,yc,srs,see,
  43.     sum_y,sum_y2    : real;
  44.  
  45. begin        { procedure linfit }
  46.   ncol:=3;    { number of terms }
  47.   for i:=1 to nrow do
  48.     begin    { setup x matrix }
  49.     xi:=x[i];
  50.     xmatr[i,1]:=1.0;    { first column }
  51.     xmatr[i,2]:=1.0/xi;        { second column }
  52.     xmatr[i,3]:=ln(xi)    { third column }
  53.     end;
  54.   square(xmatr,y,a,g,nrow,ncol);
  55.   gaussj(a,g,coef,ncol,error);
  56.   sum_y:=0.0;
  57.   sum_y2:=0.0;
  58.   srs:=0.0;
  59.   for i:=1 to nrow do
  60.     begin
  61.       yi:=y[i];
  62.       yc:=0.0;
  63.       for j:=1 to ncol do
  64.     yc:=yc+coef[j]*xmatr[i,j];
  65.       y_calc[i]:=yc;
  66.       resid[i]:=yc-yi;
  67.       srs:=srs+sqr(resid[i]);
  68.       sum_y:=sum_y+yi;
  69.       sum_y2:=sum_y2+yi*yi
  70.     end;
  71.   correl_coef:=sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow));
  72.   if nrow=ncol then nm:=1
  73.   else nm:=nrow-ncol;
  74.   see:=sqrt(srs/nm);
  75.   for i:=1 to ncol do        { errors on solution }
  76.     sig[i]:=see*sqrt(a[i,i])
  77. end;        { LINFIT }
  78.  
  79.